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Abstract 

We propose a spectral Particle-In-Cell (PIC) algorithm that is based on the combination of a Hankel transform and 
a Fourier transform. For physical problems that have close-to-cylindrical symmetry, this algorithm can be much 
faster than full 3D PIC algorithms. In addition, unlike standard finite-difference PIC codes, the proposed algorithm 
is free of spurious numerical dispersion, in vacuum. This algorithm is benchmarked in several situations that are of 
interest for laser-plasma interactions. These benchmarks show that it avoids a number of numerical artifacts, that 
would otherwise affect the physics in a standard PIC algorithm - including the zero-order numerical Cherenkov 
effect. 
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Introduction 

Particle-In-Cell (PIC) algorithms [1, 2] are extensively used in several areas of physics, including the study of 
astrophysical plasmas, fusion plasmas, laser-plasma interactions and accelerator physics. Yet, despite their wide 
use, PIC algorithms can be very computationally demanding, especially in three-dimensions, and are still subject 
to a range of numerical artifacts. These shortcomings can be particularly significant when simulating accelerated 
particle beams, or laser-plasma interactions (such as laser-wakeheld acceleration) for two reasons: 

• these systems often have close-to-cylindrical symmetry (e.g. particle beams and laser pulses are often 
cylindrically symmetric). This prevents the use of 2D Cartesian or 2D cylindrical PIC algorithms (which 
are only well-suited for slab-like and azimuthal symmetry), and is instead often dealt with by using 3D 
Cartesian PIC algorithms, which can be very computationally expensive; 

• the physical objects of interest (e.g. the laser, or the accelerated particle beam) often propagate close to 
the speed of light. This makes them very sensitive to spurious numerical dispersion , i.e. the fact that the 
electromagnetic waves do not propagate exactly at the physical speed of light in a standard PIC code, but 
travel instead at a spuriously-altered, resolution-dependent velocity. In the above-mentioned cases, spurious 
numerical dispersion can lead to substantial numerical artifacts which can mask or disrupt the physics at 
stake in the simulation. This includes, for instance, numerical Cherenkov effects in general [3], but also 
more specific artifacts, such as e.g. the erroneous prediction of the dephasing length in laser-wakefield 
acceleration [4], 

Yet several modifications can be made to the PIC algorithm, in order to mitigate these difficulties and increase the 
speed and accuracy of the simulations in these physical situations: 

• one of these modifications is the development of cylindrical PIC algorithms with azimuthal Fourier decom¬ 
position [5, 6, 7] of the electromagnetic field components (sometimes referred to as quasi-JD algorithms, 
or as quasi-cylindrical algorithms as we do here). By taking into account the symmetry of the system, these 
algorithms can typically reduce the cost of the simulation to a few times that of a 2D Cartesian simulation, 
instead of that of a full 3D Cartesian simulation. Moreover, unlike 2D Cartesian algorithms, these algo¬ 
rithms are well adapted to close-to-cylindrical physical systems and can accurately capture physical effects 
that are intrinsically 3D (such as e.g. the non-linear self-focusing of an intense laser in a plasma [8]); 
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• a second, separate modification was introduced by the development of spectral Cartesian PIC algorithms 
[9, 10, 11, 12, 13] i.e. algorithms that solve the Maxwell equations in Fourier space. (These algorithm are 
also sometimes referred to as pseudo-spectral algorithms when they make use of an intermediate interpo¬ 
lation grid in real space, and this term then also applies to the algorithm presented here.) These algorithms 
contrast with finite-difference algorithms, which solve the Maxwell equations by approximating the deriva¬ 
tives as finite differences on a discrete spatial grid. Importantly, while finite-difference algorithms suffer 
from spurious numerical dispersion, there is a class of spectral algorithms - often referred to as Pseudo- 
Spectral Analytical Time Domain (PSATD) algorithms [10, 11] - which exhibits no spurious numerical 
dispersion in vacuum. As a consequence, these algorithms are free from the associated numerical artifacts. 
Moreover, it was shown recently [14, 15, 16, 17] that spectral PIC algorithms have better stability properties 
when performing PIC simulations in a Lorentz-boosted frame [18, 19, 20]. These stability properties are 
very promising, since boosted-frame simulations can be faster than their laboratory-frame counterparts by 
several orders of magnitude [18]. 

Even though these two improvements are both very valuable, they cannot be combined with each other in 
their present formulation. Fundamentally, this is because the quasi-cylindrical algorithm uses a cylindrical system 
of coordinates, whereas most spectral algorithms (and in particular the PSATD algorithms) were developed in a 
Cartesian system of coordinates. This incompatibility is however not definitive, and the aim of this article is to 
overcome these differences by developing a spectral quasi-cylindrical formalism. In this document, we derive this 
formalism, and we show that it can be used to build a PIC algorithm that combines the speed of quasi-cylindrical 
algorithms with the accuracy and lack of spurious numerical dispersion of the PSATD algorithms. 

Although the algorithm described here is, to our knowledge, unique in its capabilities, there are several other 
existing codes which have some similarities with it. One such example is the hybrid, quasi-cylindrical version 
of OSIRIS [21], This algorithm involves a Fourier transform in the longitudinal direction, but retains a finite- 
difference formulation in the transverse (radial) direction. As a consequence, this algorithm is not fully spectral, 
nor fully dispersion-free. Another such example is the code PlaRes [22], which is designed to simulate the physics 
of free-electron lasers (FEE). This code also uses a spectral quasi-cylindrical formalism, but it models the fields 
in narrow spectral intervals around the resonant FEE frequencies, relies on the scalar and vector potentials (f> and 
A and uses no spatial grid. As a result, PlaRes differs significantly from the standard PIC formulation, and from 
the algorithm described here. 

In the present article, we start by deriving the equations of our spectral quasi-cylindrical formalism (Section 1). 
We then explain how these equations were discretized and implemented in a fully working PIC code (Section 2), 
and we report on the results of a number of benchmarks that were performed with this code (Section 3). Finally, 
we describe two typical physical situations in which our spectral algorithm performs better than standard finite- 
difference algorithms (Section 4). 

1. Representation of the fields and continuous equations 

1.1. A reminder on spectral Cartesian codes 

It is well-known that the Maxwell equations in Cartesian coordinates 

4 d t E x = d x B z - d-B x - p 0 jx d,B x = -d y E- + d z E x 

^_d,Ey = d : B x - d x B z - p {) j y d,B y = -d z E x + d x E z 

4 dt E z = d x B x - d y B x - pojz d,B z = -d x E y + d y E x 

can be solved by representing the fields as a sum of Fourier modes. 
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where F is any of the fields E, B or j, and where u is either x, y or z. T represents the Fourier components of 
F, which will be denoted S, & or J~ depending on whether F represents E, B or j. With this representation, the 
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different Fourier modes decouple and the equations Eqs. (la) to (lc) become 
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The Fourier coefficients & and B can be then integrated in time, and transformed back into real space using Eq. (2). 
This is the core principle of spectral Cartesian algorithms, including the PSATD algorithms. 


1.2. Spectral quasi-cylindrical representation 

The Fourier representation Eq. (2) is no longer the appropriate representation when the Maxwell equations are 
written in cylindrical coordinates. 
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When replacing the representation Eq. (2) into the Eqs. (5a) to (5c), the Fourier modes do not decouple. Instead 
one has to use the Fourier-Hankel representation: 


F z (r) = 
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where F is either E, B or j, where J m denotes the Bessel function of order m, and where 'F zm , r F+, m and F Jn 
represent the spectral components of F. (See Appendix A for a derivation of the above equations.) Conversely, 
the spectral components F zjn , T f .„, and F x „, are related to the real-space fields F r , F g and F z by: 
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In the above equations, notice here that the Cartesian component F z and cylindrical components F r , Fg do not 
transform in the same manner, which is due to their different behavior close to the axis. (Again, see Appendix A 
for a more detailed explanation.) Note also that scalar fields (like p) transform in the same way as the Cartesian 
component F-. 

When replacing Eqs. (6a) to (6c) into the Maxwell equations in cylindrical coordinates Eqs. (5a) to (5c), the 
different modes decouple, and the equations for the spectral coefficients become: 
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(See Appendix B for a derivation of these equations.) Notice that these equations have a similar structure as the 
spectral Cartesian equations Eqs. (4a) to (4c), as they also involve the product of the fields with the components 
of k , in their right-hand side. They do differ however in their details, as evidenced by the signs, the factors 1/2 and 
the presence or absence of the complex number i in some terms. 

Similarly, in this formalism, the conservation equations V • E - pie o and V • B = 0 become 

k ± (8+, m ~ £-.m) + ik z S z ,m = — kj_(S +tm -£-„) + ik z £ z , m = 0 (9) 

£o 

As expected, the above conservation equations Eq. (9) are preserved by the Maxwell equations Eqs. (8a) to (8c), 
provided that the current satisfies d t p + V • j = 0, i.e. in spectral space: 

&lPm "t" +j n — SF-,m) A;// . m — 0 (10) 

(This can be checked, for instance, by differentiating Eq. (9) in time and by using Eqs. (8a) to (8c) to re-express 
the time derivatives.) 

As in a spectral Cartesian PIC codes, the equations Eqs. (8a) to (8c) can be integrated in time, and the fields 
can then be transformed back into real space, using Eqs. (6a) to (6c). Therefore, the methods used to integrate the 
fields in time in a spectral Cartesian code (e.g. PSATD) should be transposable to this spectral quasi-cylindrical 
formalism. 

However, the advantage of this formalism over its Cartesian counterpart is that the sum over m in Eqs. (6a) 
to (6c) can generally be truncated to only a few terms, for physical situations that have close-to-cylindrical sym¬ 
metry. (A similar truncation is done in [6].) This is because the different values of in correspond to different 
azimuthal modes of the form and because these modes are typically zero for large values of |m|, in situations 
with close-to-cylindrical symmetry. As a result, the representation of the fields is reduced to a few 2D arrays 
T m (k z ,kJ instead of the Cartesian 3D arrays T(k x , k y , k z ). This reduction makes the manipulation of the fields 
much more computationally efficient. 


2. Numerical implementation 


2.1. Overview of the algorithm 

The PIC algorithm described here uses the above-mentioned representation of the fields, in order to solve the 
Maxwell equations and the motion of charged particles in a finite-size simulation box. 

Note that, when using spectral algorithms in a finite box, it is often necessary to arbitrarily adopt specific 
boundary conditions (which may not match the physics at stake), in order to be able to conveniently represent the 
fields. For instance, in Cartesian spectral codes, periodic boundaries are arbitrarily chosen in order to be able to 
represent the fields as a discrete sum of Fourier modes. However, this is mostly for mathematical convenience 
and it does not preclude the application, at each timestep, of another type of boundary condition in real space 
and within the finite box (e.g. Perfectly Matched Layers), before transforming the fields to spectral space (see 
e.g. [23, 24]). In the same spirit, here we arbitrarily impose periodic boundary conditions along z, and Dirichlet 
boundary conditions along r (more specifically E(r max ) = 0 and B(r max ) = 0), since in this case the fields can 
be decomposed into a discrete Fourier-Bessel series (see e.g. [25]). Yet again, this does not prevent the practical 
application of other boundary conditions in real space. 

Although, as explained in Section 1.2, the fields are represented by a few 2D arrays, the particles are still 
distributed in 3D and their motion is integrated in 3D Cartesian coordinates. As in a spectral Cartesian code, we 
do not perform the current deposition and field gathering directly from the macroparticles to the spectral space. 
(This is inefficient since the deposition and gathering are local operations in real space - i.e. affect only the few 
cells next to the macroparticle - but global operations in spectral space - i.e. they affect all the spectral modes 
simultaneously.) Instead we use an intermediate grid where these operations can be performed locally, and where 
the fields F um are defined by 
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where F is either E, B or j and u is either z, r or 9. Notice that this representation is the same as that of [6, 7]. In 
this representation, the spectral decomposition in the azimuthal direction is preserved, since the factors e"” 0 can 
be efficiently computed from the particle Cartesian positions x,y,z by using the relation e" rB = (x + iy) m /r m [6]. 
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After the particles deposit their charge and current onto this intermediate grid, the fields of the grid are trans¬ 
formed into spectral space where, as mentioned in Section 1, the Maxwell equations can be easily integrated. 
From Eqs. (11) and (12) and Eqs. (6a) to (6c) and (7a) to (7c), the transformation between the intermediate grid 
F m (r,z) and the spectral grid 'f r m (k ± , k z ) is: 


and 


t z Jk x ,k z ) = HT,„[ FT[ F z , m (r, Z ) ] ] 
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T- m (k ± ,k z ) = HT m _i 
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F z .Jr,z) = IFT[ IHT m [f z , m (k ± ,k z )] ] 

F rjn (r,z) = TFT| lHT m+l [f + . m (k ± ,k z )] + IHT m _i [f'- m (k ± ,k z )] ] 
F g ,m(r, z) = i TFT [ im m+l [f + . m (k ± ,k z )] - [f-, m (k ± ,k z )] ] 


(13a) 
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(14b) 

(14c) 


where FT represents a Fourier Transform along the z axis and HT„ represents a Hankel Transform of order n along 
the transverse r axis, and where IFT and IHT„ represent the corresponding inverse transformations: 
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FT [f](k z )= d ze~ ik *f(z) IFT[g] (z) = — d k z e^~ g(k z ) (15) 
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The above equations show that the transformation from the intermediate grid (F u m ) to the spectral grid ( r F u , m ) 
is the combination of a Fourier transform (in z) and a Hankel transform (in r). The Fourier transform in z can 
be discretized through a Fast Fourier Transform (FFT) algorithm, which requires an evenly-spaced grid in z and 
in k z . On the other hand, there is more freedom of choice for the Discrete Hankel Transform (DHT), and the 
implementation that we chose is described in the next section (Section 2.2) and in the Appendix D. 

Figure 1 gives an overview of the successive steps involved in one PIC cycle, including the respective role of 
the intermediate grid (F um ) and spectral grid ('F lL ,„). Note that, in the PSATD scheme that we chose (and which 
is described in more details in Section 2.6), all the fields are defined at integer timesteps, except for the currents, 
which are defined at half timesteps. Sections 2.3 to 2.6 describe the successive steps of the PIC cycle in more 
details. 


2.2. Transverse discretization of the intermediate grid, of the spectral grid and of the Hankel Transform 

When transforming from the intermediate to the spectral grid, the FFT algorithm is the natural way to discretize 
the Fourier transform along z, due to its favorable computational scaling (oc N z log (N z )). On the other hand, there 
are a variety of existing algorithms (that are not mathematically equivalent) to discretize the Hankel transform 
(e.g. [26, 27, 28, 25, 29]), and the use of one or the other is very dependent on the application pursued. Broadly 
speaking, choosing an algorithm for the Discrete Hankel Transform (DHT) algorithms consists in: 

• choosing a discrete grid in r and k ± space, on which to sample the functions to be transformed. In some 
algorithms, these grids may not be evenly-spaced, and can be for instance logarithmically spaced [28] or 
can correspond to the zeros of Bessel functions [27, 25, 29], 

• once a grid is chosen, the DHT amounts to a linear operation on a finite set of points (the points of the 
grid) and can thus be represented by a matrix. Thus the second choice is that of a matrix that represents, as 
closely as possible, the exact Hankel Transform. 

Here, we choose to discretize the algorithm on an evenly-spaced grid in r (for the intermediate grid) : 

= ; e [0,.... N r - 1} where A r = ^ (17) 

but on an irregular grid in k± (for the spectral grid). This is because, for the chosen boundary condition (E(r max ) = 
0, B(r max ) = 0), the fields can be expressed as a discrete sum of Bessel modes J m with 


j e {0. N r - 1} 
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b) Equations of motion 
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c) Current and charge deposition Current correction 



d) Maxwell equations 
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Figure 1: Schematic description of the 4 steps of a PIC cycle. At any given time, the quantities that are known 
are shown in color (red and blue), while the quantities that are unknown or have been erased from memory are 
shown in gray. The quantities that are being calculated at a given step are displayed with a colored background, 
and arrows indicate which quantities are used for this calculation. 
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where a” 1 is the jth positive zero of the Bessel function of order m J m (including the trivial value a"J = 0 for m > 0; 
see e.g. [25]). These values are represented in figure 2. Notice that it is not an issue that the spectral components 
r F m for different azimuthal modes m are discretized on different k ± grids, since each azimuthal mode m evolves 
separately in Eqs. (8a) to (8c). 



Index j Index j Index j 


Figure 2: Position of the grid points in k± space (blue dots) for the azimuthal modes in = 0, m = 1 and m = 2 for 
N r = 10 and r max = 1. These values are compared with those of an evenly-spaced grid used in spectral Cartesian 
codes (grey dots, k, = jn/r max ). 

As mentioned above, once this grid is set up, the Discrete Hankel Transform is simply a linear operation on a 
finite set of points, and can thus be represented by a matrix operation. 


JV r -l JV r -l 

DHT”[/] (k’lj) - 2 (M nm )j p f(r p ) IDHT”[g] (r ; ) - j] (M' lum )j, P g(k%) (19) 

p =0 p =0 

Notice that the N r x N r transformation matrices M njn and M' n m depend on n (order of the Hankel transform, 
i.e. order of the Bessel function J„ in the integrand of Eq. (16)) and in (index of the azimuthal mode, and thus 
index of the spectral grid k r " ■ on which the Hankel transform is performed). In practice, n and m are either equal 
or they differ by +1 (e.g. in Eqs. (14b) and (14c), when the azimuthal mode in is transformed using the Hankel 
transform of order m + 1 or m - 1.) 

The expressions of M n m and M' n m , for the DHT algorithm that we chose, are given in Appendix D. In practice, 
these matrices need to be computed only once (at the beginning of the simulation) and can then be used at each 
iteration. Notice also that the computational time for this matrix multiplication is proportional to N which is 
slow compared to a ID FFT (oc N r \og(N r )), but still faster than the 2D FFT (oc N x N y \og{N x ) + N x N y \og(N y )) 
which is typically used in the transverse plane of a spectral 3D Cartesian code. 


2.3. Field gathering 

In the following, we describe the four successive steps of a PIC cycle with our algorithm in detail, starting 
with the field gathering. When gathering the fields from the intermediate grid to the macroparticles (see Fig. 1), 
we use the standard linear shape factors : 
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( 20 ) 


where F is either E or /(, u is either r, 8 or z, k is the index of the macroparticle, and p and q are the indices of 
the two nearest cells in z and r respectively. N m is the total number of azimuthal modes used, and % denotes the 
real part of a complex quantity. Notice that, in Eq. (20), we used the fact that F u _ m {z, r) = F* m (z , r), which can be 
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inferred from Eq. (12). Incidentally, Eq. (20) shows that only the modes with m > 0 need to be taken into account 
in the code, as they are sufficient to retrieve the force on the macroparticles. 

Finally, in Eq. (20), S ,,, and S are the linear shape factors in z and r : 


Zp+ 1 z 

A z 

S z , p+ i(z) = 

l <1 

with z p < z < z p +i 

(21a) 

r q + 1 - r 

S r . q +i(r) = 

r-r q 

with r q < r < r q+l 

(21b) 

A r 

A r 


For particles that are in the lower half of the first radial cell (r < r o = Ar/2), Eqs. (20) and (21b) require the value 
Fu,m(Zp , r \) although r \ = -Ar/2 is not actually part of the grid. In order to still apply Eqs. (20) and (21b), we 
explicitly set F um (z p ,-Ar/2) = ±F ujn (z p , Ar/2), where the - sign is chosen whenever the considered field F up , is 
by definition zero on the axis, and the + sign is chosen otherwise, (e.g. F r ,a is by definition zero on the axis ; see 
[6] for more details, and for a similar method.) 

Note that we use linear shape factors here only for the sake of simplicity, and that higher-order shape factors 
(e.g. quadratic or cubic) could also be used in principle, with a similar mirroring of the field values across the axis. 


2.4. Equations of motion 

Since the macroparticles evolve in 3D, we first compute the Cartesian components E x , E y , E z , B x , B y and 
B z , from the fields E r , Eg, E z , B r , B„ and B that were gathered at the positions of each macroparticle. We then 
advance the equations of motion 


dp dx p 

— - qE + qv X B — — - 

dt dt ym 

in standard 3D Cartesian coordinates by using the leap-frog pusher described in [30], 


( 22 ) 


2.5. Current deposition 

As in [6], the charge density is calculated on the intermediate grid in the following way: 

, Z k S z , p (z k )S r , q (r k )Q k e^ 

Pm\Zp, fq) — 

*q 


(23) 


where Qj ; is the charge of the macroparticle with index k, and where V q is the volume of a cell, which is for our 
grid 

V q = 7t[(q+l) 2 -q 2 ]Ar 2 Az (24) 

Similarly, the current deposition is given by 


, , , Tik S Z ,p(Zk)S r,q(j'k)Q k Vu,ke imek 

Ju,m(z p , r q ) = -—- (25) 

v q 

where u = z,r, 8 and the v u k are the cylindrical components of the velocity of the macroparticle k. As mentioned 
previously, once the macroparticles have deposited their charge and current on the intermediate grid, we transform 
them to the spectral grid, using an FFT and a DHT (see Fig. 1). 

Note that, when a macroparticle travels through the grid, the factor V q in the above formulas decreases as the 
macroparticle comes closer to the axis. Therefore, close to the axis, the impact of a single macroparticle on an 
individual grid cell can be substantial and this generally translates into higher levels of noise. (This is a general 
problem with cylindrical and quasi-cylindrical codes, and thus it is not specific to the present spectral version.) In 
order to mitigate this problem, and more generally in order to avoid the accumulation of noise at high frequency, 
smoothing is typically applied on the charge and currents, after they have been deposited (but no smoothing is 
applied on the fields E and B directly). This smoothing is performed directly in spectral space, by multiplying the 
charge and currents by a transfer function T(k : , k ± ) which damps the high frequencies k. and whose mathematical 
form is identical to the spectral Cartesian representation of a single-pass binomial filter [1], 

t(k z ,k r ) = c °s 2 ( T~~ cos 2 ( t ~— (26) 

\ &z,max Z ) \ K±,max Z ) 

where k- ynnx and k ±tinax are the highest wavevectors that the discrete spectral grid supports, in the longitudinal and 
transverse direction. 

Notice also that, in the above deposition scheme, we do not attempt to reproduce the Esirkepov charge- 
conserving current deposition [31], and instead use a simple direct current deposition. It is well-known that 
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this simple deposition does not necessarily satisfy the relation d t p + V • j = 0, but that this can be corrected, by 
slightly modifying the currents without modifying their curl (e.g. [32]): 

/ = j ~ VG (27) 


where G satisfies the Poisson-like equation 

V 2 G = d t p + V • j 


(28) 


The above equation is typically expensive to solve on a spatial grid, but very easy to solve in spectral space. 
In spectral space and with the notations of Fig. 1, these equations become 


fy-t n+l/2 A-n+\/2. k ±. An+l/2 

+,m +,m ' 2 ^ m 


fj-rn+1/2 _ Sr-n+l/2 ^-L ^A/i+1/2 

-,m ~ -,m 2 ^ m 


rrf n+ 1/2 _ rj-n+l/2 _ •/ An+l/2 

'J z,m — z,m ll ^z&m 


(29) 


with 
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At 
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(30) 


With this correction, the new currents J~ ,n+1 / 2 do satisfy the charge conservation equation Eq. (10). Therefore we 
apply this correction in spectral space, at the end of the current deposition at each timestep. 


2.6. Integration of the Maxwell equation using the PSATD scheme 

The Maxwell equations Eqs. (8a) to (8c) could in principle be integrated by using a finite-difference scheme 
in time, which would be an adaptation of the Cartesian PSTD scheme [13], 
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(32a) 

(32b) 

(32c) 

(33a) 

(33b) 

(33c) 


However, this type of scheme retains some amount of spurious numerical dispersion, and can thus affect the 
simulated physics. 

Instead, here we use an adaptation of the Cartesian PSATD scheme [10] for our spectral quasi-cylindrical 
representation. As in the case of the standard PSATD, we assume that the currents are constant over one timestep 
and that the charge density is linear in time over the same timestep. Under these assumptions, the Maxwell 
equations Eqs. (8a) to (8c) can be integrated analytically over that timestep, and they lead to: 



(34a) 

(34b) 

(34c) 
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(35c) 


- (M'; m + ^£", m )+/i 0 c 2 —— (ik.jiT + ik^jiT) 

Cl) ' 7 Cl) z x ' 

where w = c yk? + C = cos(wAf) and 5 = sin(ruAr). (See Appendix C for a derivation of these equations.) 


2.7. Practical implementation 

The full PIC algorithm described in this section was implemented in the code FBPIC (Fourier-Bessel Particle- 
In-Cell), which is written in Python. For performance, this implementation makes use of the pre-compiled libraries 
FFTW [33] and BLAS [34] (for the matrix multiplication in the DHT), and utilizes the Numba just-in-time com¬ 
piler [35] for the computationally-intensive parts of the code (current deposition and field gathering). In addition, 
the code was developed for both single-CPU and single-GPU architectures (with the GPU runs being typically 
more than 40 times faster than the equivalent CPU runs, on modern hardware). Importantly, a precise timing of 
the different routines showed that, although the time taken by the spectral transforms (FFT and DHT) is not entirely 
negligible, it does not usually dominate the PIC cycle. For instance, on a K20 GPU and for a grid with N z = 4096, 
N r = 256 and 16 particles per cell, the FFTs and DHTs take up 6% and 14% of the PIC cycle respectively, while 
the rest of the time is dominated by the field gathering and current deposition. 

Notice that, even though parallelization is generally challenging for spectral algorithms, a multi-CPU/multi- 
GPU version could still be developped in the future by using the method of [32], For the present algorithm, 
this would involve a domain decomposition along the z axis, whereby FFTs would be performed locally within 
each subdomain and whereby a large number of guard cells would be used in order to mitigate the errors at the 
border between subdomains. This method based on local FFTs has been studied in the Cartesian context [23], and 
is now rather well understood. By contrast, performing domain decomposition in the r direction would be more 
challenging, due to the absence of previous work on local Hankel transforms. At any rate, the above considerations 
for parallel implementation are out of the scope of the present article, and will be the subject of future work. The 
benchmarks described in the following section use a single-CPU/single-GPU version of the algorithm. 


3. Benchmarks 

We tested the above algorithm in a number of physical situations. These tests included standard problems, 
such as e.g. a periodic plasma wave, and involved comparing simulated results to analytical solution, as well as 
verifying that the algorithm conserves the energy to a satisfying level. For the sake of conciseness, in the present 
article, we restrict the discussion to the tests of the dispersion relation, since the absence of spurious numerical 
dispersion was one of the original goals of this algorithm. 

3.1. Propagation in vacuum 

A first test consisted in letting a laser pulse propagate in vacuum and in measuring its group velocity in the 
simulation. These simulations were run with the PSATD quasi-cylindrical algorithm presented in Section 2 (see 
esp. Eqs. (34) and (35)), but also, for comparison, with a PSTD version of this same algorithm (see Eqs. (31) 
to (33)), as well as with a finite-difference quasi-cylindrical algorithm equivalent to that of [6, 7] and which had 
been previously implemented in the PIC code Warp [36]. 

The simulations were run with a moving window, in a box with a longitudinal size of 40 pm and a transverse 
size of 48 pm. (In the case of the spectral algorithm, the fields were damped at the back of the moving window, in 
a similar way as in [37], in order to prevent periodic wrapping of the fields). The laser pulse itself was initialized 
at focus, with a waist wq = 16 pm, a length L - 10 pm, a wavelength A = 0.8 pm and a dimensionless potential 
vector an = 10 2 . The resolution was varied while keeping the same cell aspect ratio (Ar = 5As). The timestep 
was set to cAr = As in the case of the PSATD algorithm, to cAt = 0.9 x 2/n^l/Az 2 + 1/A r 2 in the case of the 
PSTD algorithm, and to cAt =1/ V1 /As 2 + 2/Ar 2 in the case of the finite-difference algorithm (for the PSTD and 
finite-difference algorithms, the chosen Ar is close to the Courant limit). 

Physically, the on-axis group velocity of the pulse should be slightly lower than c due to the finite waist of the 
pulse (see e.g. [38]). The analytical expression for the corresponding relative difference in on-axis group velocity 
is 

c - v e I A 

= 2 (^ 

10 


C 


(36) 










In the case at hand {wq = 16 /mi, A = 0.8 pm), this relative difference is extremely small ((c-v g )/c = 1.27x 10 4 ), 
and it may be difficult for PIC codes to capture this small physical difference. 

To assess the capacities of the finite-difference and spectral codes in this regard. Fig. 3 displays the relative 
difference in group velocity, as measured in the simulations. As can be observed, in the finite-difference simula¬ 
tions and PSTD simulations, the group velocity of the laser depends on the resolution, due to spurious numerical 
dispersion. Moreover, in these cases, the group velocity is considerably different than the analytical prediction, 
even for a relatively high resolution. (Since Fig. 3 does not display the sign of c—v g , it is worth noting here that the 
PSTD algorithm leads to v g > c while the finite-difference algorithm has v g < c. It is also important to realize here 
that, although the finite-difference algorithm performs better than the PSTD algorithm in this particular example, 
this is not generalizable to other cases, as e.g. the numerical dispersion of both algorithms behave differently when 
At or the ratio Az/Ar is changed.) 



Figure 3: Relative difference between c and the group velocity of a laser pulse in vacuum, for different resolutions. 
The dashed black line represents the analytical prediction given by Eq. (36), while the blue and red points represent 
the results of PIC simulation, with the PS ATD, PSTD and finite-difference algorithm. 


On the other hand, with the PSATD algorithm, the group velocity is practically independent of the resolution, 
and displays very good agreement with the analytical prediction. This corroborates the fact that the PSATD quasi- 
cylindrical algorithm described here has no spurious numerical dispersion in vacuum. 

Since the PSTD algorithm is less acurate than the PSATD algorithm, while not providing any substantial 
advantage in terms of speed or practical implementation, we do not consider it further here and perform the rest 
of the tests with the PSATD algorithm. 


3.2. Linear propagation in a plasma 


In order to confirm that the spectral algorithm also performs well in the presence of a plasma, we ran the same 
type of simulations with a uniform, pre-ionized plasma. The numerical parameters of the simulations, as well as 
the physical parameters of the laser, were the same as in the previous subsection (Section 3.1). The plasma was 
represented by 16 macroparticles per cell, and had a density n e = 1.75 x 10 18 cm -3 = 10 3 n ci where n c is the 
critical density for A = 0.8 pm. 

Since the intensity of the laser is very low here (ao = 10 2 ), the propagation is linear, and the on-axis group 
velocity is given by (e.g. [38]) 


c~Vg _ 
c 2 n c 



(37) 


In Fig. 4, we compare this analytical prediction with the group velocity in the finite-difference and spectral 
simulations. Again, the group velocity is resolution-dependent in the finite-difference algorithm, and it is sub¬ 
stantially different than the analytical prediction even at high resolution. In the spectral code, the group velocity 
velocity exhibits a very weak dependence on resolution (which is likely due to the errors of current deposition and 
field gathering on the finite grid). However, its value remains always very close to the analytical prediction at all 
resolutions. 

We emphasize that, although the difference between c and v g is small here and may thus seem unimportant, it 
is this difference which determines the dephasing length and thus the maximum beam energy, in a laser-wakefield 
simulation. It is therefore paramount to obtain its correct value in the simulation codes. This problem is well- 
known in the case of standard finite-difference codes. For Cartesian finite-difference codes, a common solution 
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Number of gridpoints per laser wavelength (along z) 


Figure 4: Relative difference between c and the group velocity of a laser pulse in a plasma at 10 3 n c , for different 
resolutions. The dashed line represents the analytical prediction given by Eq. (37), while the blue and red points 
represent the results of PIC simulations. 


is to use a scheme which is dispersion-free along the z axis (e.g. [39, 40, 41]), but no such scheme has been 
developed for quasi-cylindrical codes. Alternatively, spurious numerical dispersion is often dealt with in finite- 
difference codes by either using an even finer grid in z (which is very computationally expensive) or a coarser 
grid in r. (Recall that, in the simulations shown here, A r = 5A z. A coarser resolution in r allows to use a 
slightly larger timestep At that approaches cAz, the limit at which spurious numerical dispersion vanishes in the 
^-direction.) However, a coarser resolution in r may not always be adapted to resolve the physics of interest, while 
a finer resolution in z is expensive. It is therefore remarkable that, with the spectral quasi-cylindrical algorithm, 
the correct group velocity is obtained independently of the cell aspect ratio, without having to specifically adapt 
the resolution in r or z. 


3.3. Linear laser-wakefield 

In order to further ascertain that the spectral algorithm described here gives appropriate results beyond simple 
dispersion tests, we ran a simulation of a laser-wakefield, and compared the amplitude of the wakefield with the 
corresponding analytical predictions. 

In these simulations, the laser was linearly polarized along the transverse x direction, and had an amplitude 
ao = 10~ 2 , a wavelength A = 0.8 /mi, a length L = 10 /J.m and a waist wq = 20 //m. The plasma was preionized, 
with a uniform electron density n e = 1.75 x 10 18 cm~ 3 = 10 ' n c . As is often the case for simulations of intense 
laser-plasma interaction (where thermal effects are generally assumed to be negligible), the initial temperature of 
the plasma is set to 0. The simulation was run in a moving window whose longitudinal and transverse sizes were 
80 /on and 60 /on. The spatial and temporal resolutions were A z — 0.05 /mi, A r = 0.5 /mi and At = Az/c. 
Since here a {) — 10 2 , the analytical wakefield is given by the linear and quasistatic theory. For a laser of the form 
a{z, r, t) = adz, t)e~ r ‘ w o cos(koz - «o0 - where adz, t) is the longitudinal envelope of the laser - the longitudinal 
and transverse fields E z and E y are given by (e.g. [8]) 


E z (z, r, t ) = 


r 2 kl 


r 


E,U.r. ,) = -^ k 4 
e wz 


I 


a 2 dz',t)e r2/tl »cos [k p (z-z!)]dz' 

oo 

a 2 (z\ sin[k p (z - z)]dz 


(38a) 

(38b) 


where k p is the plasma wavevector. 

Here, the longitudinal laser envelope ai was extracted directly from the simulation, and the above analytical 
integrals where carried out numerically. The resulting predicted fields E z and E y are plotted in dashed lines in 
the left and right lower panels respectively of Fig. 5. These predicted curves are compared with the fields E z and 
E y extracted directly from the simulation (red lines in the lower panels). These fields from the simulation are 
also displayed as colormaps in the upper panels of Fig. 5, and these colormaps include the line along which the 
quantities in the lower panels are plotted (dashed line). 

As can be seen on the lower panels of Fig. 5, the analytical and simulated curves overlap very precisely, which 
confirms the validity of our PIC algorithm. This agreement is all the more remarkable as the E z field has been 
extracted from the cells that are closest to the axis, a region where quasi-cylindrical algorithms are typically very 
noisy. 
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Figure 5: Upper panels: Colormaps of the fields E z and E v as extracted from the simulation (blue and red), along 
with the laser envelope (orange contour lines). (The laser propagates to the right and is polarized along x.) The 
dashed lines indicate the position where the curves of the lower panels have been extracted. Lower panels: profile 
of the fields at a given radial position, as given by Eqs. (38a) and (38b) (dashed line) or as given by the simulation 
(red line). 


4. Advantages over finite-difference algorithms 

Since finite-difference algorithms and spectral algorithms are considerably different, each of them has its own 
advantages and shortcomings. 

For instance, spectral algorithms - including the one described here - are generally more difficult to parallelize 
than finite-difference codes. Similarly, the implementation of specific boundary conditions, and of a moving 
window, is usually more challenging in a spectral algorithm. On the other hand, due to their high accuracy, spectral 
algorithms can avoid a number of numerical artifacts that are typically present in finite-difference algorithms. This 
is quite important as, in a finite-difference simulation, these artifacts may remain unnoticed unless additional care 
is taken in their analysis, and yet they may affect the physics at stake in an important way. 

One example of such an artifact is spurious numerical dispersion, which, as mentioned in [4] and in Section 3.2, 
can modify the dephasing length of a laser-wakefield accelerator - a spurious effect which may be difficult to 
discern, especially in the cases where there is no analytical formula for this length. As shown in Section 3.2, the 
algorithm described here would not suffer from this effect, as it is free of spurious numerical dispersion. Similarly, 
in this section we describe two other typical situations in which our spectral quasi-cylindrical algorithm avoids 
important artifacts, that would otherwise arise in a finite-difference code. 

4.1. Suppression of zero-order numerical Cherenkov effect 

A first artifact which is avoided is the zero-order numerical Cherenkov effect. To zero order, the numerical 
Cherenkov effect is a consequence of the spurious numerical dispersion [3], and arises because some relativistic 
particles can travel faster than the numerically-altered velocity of the electromagnetic waves, in the simulations. 
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This typically causes these relativistic particles to emit a characteristic spurious radiation. Recently, it was shown 
that this spurious radiation can have a very substantial impact in lab-frame simulations of laser-wakefield accel¬ 
eration, in particular by spuriously increasing the emittance of the accelerated beam [42], It was also shown that, 
by modifying the PIC algorithm and its numerical dispersion relation, the zero-order Cherenkov radiation can be 
suppressed. (Notice that there exists also a set of higher-order, aliased numerical Cherenkov effects, which are 
not as easily suppressed [43, 44, 14, 15, 16, 17, 45]. These high-order effects can lead to disruptive instabilities in 
boosted-frame simulations, but on the other hand no such instability was observed in typical lab-frame simulations 
of laser-plasma acceleration. It is in fact likely that these instabilities may not have enough time to develop in the 
case of typical lab-frame simulations.) 

Since the spectral quasi-cylindrical algorithm described here is dispersion-free, it should not exhibit the zero- 
order numerical Cherenkov effect. In order to confirm this prediction, we ran a simulation of laser-wakefield 
acceleration with our spectral quasi-cylindrical algorithm, and compared it with an equivalent simulation that uses 
the finite-difference quasi-cylindrical algorithm of Warp. 

In these simulations, a laser pulse with a waist wo = 16 pm, a length L - 10 pm, an amplitude a () = 4 and 
a wavelength A = 0.8 pm is sent into a longitudinally-shaped pre-ionized gas jet. The gas jet has a 100 pm-long 
rising density gradient at its entrance, followed by a 200 pm plateau at a density n e = lx 10 18 cm~ 3 and then by 
a 100 pm-long downramp, so as to finally reach a density n e = 0.5 x 10 18 cm~ 3 . Again, the initial temperature 
of the plasma is set to 0. The downramp causes the injection of an electron beam, which is then accelerated in 
the subsequent density plateau at n e = 0.5 x 10 18 cm -3 . Both the spectral and the finite-difference simulations 
were run in a moving window whose longitudinal and transverse dimensions were 160 pm and 48 pm, and both 
simulations used a resolution A z = 0.032 pm and A r — 0.19 pm. Again, the finite-difference simulation was run 
with a timestep cAt - 1 / 1 / A~ 2 + 2/Ar 1 while the spectral simulation was run with cAt = Az. Importantly, the 

charge and currents are smoothed in both simulations. The spectral algorithms performs smoothing in spectral 
space as described in Section 2.5, while the finite-difference algorithm applies a single-pass binomial filter on the 
spatial grid, in both the z and r directions. 

Figure 6 shows a snapshot of the two simulations, at a similar physical time. The red colormap represents the 
quantity \E y + cB x \, while the superimposed shades of blue represent the electron density. (The quantity E y + cB x 
is chosen because it corresponds to the y component of the Lorentz force F = qE + qv x B felt by a relativistic 
electron having v = ce z .) As can be seen in the top panels, the global aspect of the bubble and of the injected 
bunch is similar in both simulations. In particular, the injected charge was found to be almost the same (766 pC 
with the finite-difference algorithm and 750 pC with the spectral algorithm). However, a closer look at the bunch 
in the finite-difference algorithm (see the middle left panel in Fig. 6) reveals that the bunch emits a high-frequency 
radiation. This radiation seems to be due to the zero-order numerical Cherenkov effect - an interpretation which 
is confirmed by the observation of the corresponding characteristic double-parabola in k space [42, 43, 44] on 
the lower left panel of Fig. 6. Importantly, the amplitude of this unphysical field, in the middle panel of Fig. 6, 
happens to be comparable to the focusing fields inside the bubble, and it can thus potentially affect the bunch in 
a substantial manner. On the other hand, the spectral algorithm does not exhibit this unphysical radiation in real 
space (see the middle right panel in Fig. 6), nor does it exhibit the corresponding characteristic pattern in k space 
(see the lower right panel). This was indeed expected from its dispersion-free property, and it confirms the fact that 
our spectral quasi-cylindrical algorithm is free of the unphysical artifacts associated with the zero-order numerical 
numerical Cherenkov effect. 

For completeness, we remark that a similar suppression of the numerical Cherenkov radiation was obtained in 
[42, 4], by using modified finite-difference algorithms. However, it is important to note that these algorithms were 
applicable to a Cartesian PIC code, but not to a quasi-cylindrical code. An adaptation to cylindrical geometry is 
proposed in [46], and was shown to efficiently suppress zero-order numerical Cherenkov effect. However, this 
adaptation is not, strictly speaking, dispersion-free. 

4.2. Accurate force of a laser on a copropagating electron 

Another case in which our spectral algorithm performs better than finite-difference algorithms is whenever a 
relativistic electron bunch overlaps with a copropagating laser. This situation occurs for instance in laser-wakefield 
acceleration, when the accelerated electron bunch progressively catches up with the laser pulse and may eventually 
overlap with it (e.g. [47, 48]). It is also the typical configuration for simulations of free-electron lasers, where 
the overlap of the electron bunch and the copropagating laser radiation, inside an undulator, leads to a growing 
instability. 

Despite the practical importance of this physical configuration, it was recently shown that standard finite- 
difference PIC codes tend to largely overestimate the force felt by the electrons inside a copropagating laser (see 
the appendix of [49]). This overestimation was shown to be mainly due to the staggering in time of the E and B 
fields in a standard finite-difference code. This staggering indeed results in an improper numerical compensation 

14 



30 


Finite-difference 


Snectral (PSATD1 




k . (pm 1 ) 


Cl 

o 


+ 


CM 

H 

Uh 

Cm 


SP 



620 640 660 680 

* (/an) 



632 


3 


634 636 

z (/tm) 


638 



.80 


72 


64 

1 

56 

48 

> 

O 

40 


32 

24 

H 

oq 

o 

+ 

16 

KL 

8 


0 


.80 


72 


64 

| 

56 

48 

> 

o 

40 


32 

24 

H 

K) 

o 

+ 

16 

kT 

8 


0 


I 15 


L 

> 

13 

oq 

u 

+ 

12 


11 

CM 

H 

Uh 

lio 

o 

i 


Figure 6: Snapshot of a finite-difference simulation (left) and of a spectral simulation (right), at a similar physical 
time. Upper panels: Representation of the bubble in real space; the superimposed blue shades represent the 
electron density. Middle panels: More detailed view of the region surrounding the bunch. Lower panels: two- 
dimensional Fourier transform of the quantity E y + cB x ; the central spot around k z = 0 corresponds to the physical 
fields, while the double-parabola in the left panel is a signature of the zero-order Cherenkov effect. 


of the two terms of the Lorentz force F = qE + qv x B. However, in our spectral quasi-cylindrical algorithm, the 
fields E and B are not staggered in time, and thus the force on a copropagating electron bunch should be correct. 

To confirm this prediction, we ran a test simulation in which a single relativistic macroparticle copropagates 
with a laser pulse. The initial configuration of the simulation is represented in Fig. 7. The laser pulse is polarized 
along x and is characterized by a waist wo = 25 pm, a length L = 7 pm, an amplitude a o = 0.2 and a wavelength 
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A = 0.8 pm, while the macroparticle represents a relativistic electron having an initial Lorentz factor y e = 25. 
For comparison, the simulation was run, again, with both our spectral quasi-cylindrical algorithm and the finite- 
difference quasi-cylindrical algorithm of Warp. In order to study numerical convergence, the simulations were 
mn with various longitudinal resolution A z, but with a fixed cell aspect ratio A r = 10Az. The timestep was again 
cAt — 1 / -\/1 / Az 2 + 2/ Ar 2 for the finite-difference simulation and cAt = A; for the spectral simulation. Finally, the 
simulations were run in a moving window with a longitudinal and transverse size of 40 pm and 60 pm respectively. 



£ (pm) 


Figure 7: Representation of the situation which is simulated in Section 4.2. The red and blue colormap corresponds 
to the electric field of the laser pulse, while the black dot corresponds to the electron considered. The laser 
propagates to the right, and is polarized along x. 

In order to assess the correctness of the force felt by the electron, we analyze the evolution of its transverse 
momentum in time. Analytically, it can be shown from the Lagrangian £, — - y 1 - f} 2 me 2 — ev ■ A that the 
equation of motion along the x axis (i.e. the axis of laser polarization) is 

1 d l p x \ p x da x a 2 0 

- —- a x ) = -—- (39) 

c at \ m e c ) ym e c ox ywo 

where a x is the normalized vector potential of the laser pulse. The right-hand side corresponds to the ponderomo- 
tive force, and, for the parameters and timescale considered here, it is in fact negligible. The canonical momentum 
of the electron is thus approximately conserved. 


Px 

—--a v = const. (40) 

m e c 

Since, inside the laser pulse, a x oscillates between the values —ao and ao, the quantity p x /m e c is predicted to also 
also oscillate between -ao and ao, as the electron progressively dephases with the laser. 

Keeping in mind that here initially an = 0.2, one can observe on Fig. 8 that the oscillation amplitude of 
p x /m e c is much higher than analytically predicted in the finite-difference simulation. In addition, these oscillations 
strongly depend on the resolution of the simulation, which confirms their spurious nature. These observations are 
consistent with those of [49], and they are due to the above-mentioned erroneous calculation of the Lorentz force. 
On the other hand, in the spectral simulations, the oscillations of p x /m e c exhibits virtually no dependence on the 
resolution. In addition, the oscillations of p x /m e c have a realistic amplitude. (The fact that these oscillations do 
not reach exactly the value 0.2 may be due to the fact that ao decreases as the laser propagates, as a consequence 
of diffraction.) These observations confirm that the spectral quasi-cylindrical algorithm properly calculates the 
Lorentz force on the electron. This is also consistent with our expectations, which were based on the fact that E 
and B fields are not staggered in time, in this algorithm. 
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Propagation distance z (/(in) Propagation distance z (//in) 

Figure 8: Evolution of the transverse normalized momentum of a macroparticle, as it copropagates with a laser 
pulse having a o = 0.2 (see Fig. 7 for a schematic representation of the situation simulated). The left and right 
plots correspond to the finite-difference and spectral algorithms respectively (note the different vertical scales on 
the two plots). The colored curves correspond to different resolutions for the simulations. 


Conclusion 

In this article, we derived the equations of a spectral quasi-cylindrical PIC algorithm, and discussed its numer¬ 
ical implementation. As explained in the text, because this algorithm is quasi-cylindrical, it requires much less 
computational time and memory than a full 3D algorithm. In addition, the fact that it is spectral (and uses the 
PSATD algorithm) prevents a certain number of artifacts, that are usually present in a finite-difference algorithm. 

This second aspect was tested in detail in this article, through several benchmarks. By comparing simulations 
with analytical results, we showed that our spectral quasi-cylindrical code is free of spurious numerical dispersion, 
that it does not exhibit zero-order numerical Cherenkov effect in typical lab-frame simulations, and that it can 
accurately calculate the force of a laser on a copropagating electron. By running the same benchmarks with a 
finite-difference quasi-cylindrical algorithm, we showed that these numerical artifacts can, on the contrary, be 
present in finite-difference PIC algorithms. 

However, it must not be forgotten that finite-difference algorithms have advantages of their own, including 
easier parallelization and application of boundary conditions. In the case of our spectral quasi-cylindrical algo¬ 
rithm, parallelization to several nodes could nonetheless be achieved by following for instance the method of [32], 
This will be the subject of future work. 
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Appendix A. Derivation of the spectral quasi-cylindrical representation 

In order to derive the representation Eqs. (6a) to (6c), we have to distinguish the Cartesian components (e.g. 
E x , E y , E z , B x , B v , B), which are well-defined everywhere in space and thus have a regular Fourier representation, 
from the cylindrical components (e.g. E r , Eg, B r , Bg), which are ill-defined at r - 0. (For instance, for a field of 
the form E = Ege x , the expression of E, is E, = Eg cost/9), which ill-defined at r - 0, as 6 itself is ill-defined at 
this position.) 


Appendix A.l. Cartesian components 

Fet F u be the Cartesian component of a field F (typically F is E, B or / and u is x, y or z). Its Fourier 
representation is thus given by Eqs. (2) and (3): 


F u (r) = 


(2/r) 3 


s-*oo s-*oo s-*oo 

Idk x Id ky Idk z T u (k)e ,( F x+k >y +k ’E> 

— oo J —OO \J —CO 


/~*oo /~*oo /~*oo 

T u (k) = dx I dv I d- F„ 

%J — OO kJ —OO — OO 


(r)e 


-i(k x x+k y y+k z z) 


(A.l) 
(A.2) 


Using the change of variable k x = k ± cos ((f)), k y = k ± sin(^), x = rcos(ff), y = r sin(0), this becomes 

1 poo /'■‘Oo p2 n 

F " (r) = (2n? J dk: J 0 k±dk± Jo d<P ^(k)e i(k±rcos(e -* )+kzZ) (A.3) 

/'-‘OO /'-‘CO /~*oo 

T u (k) = Id? j rdr Id(9 F u (r) e - i(k±rcos ^ +k ^ (A.4) 

— OO 0 sj—OO 

Fet us now use the relation e ‘Ercos(e-i/>) _ i m J m (k ± r)e"" ( ^~ (l> (which is simply another way of writing the 

well-known relation e ! “ sln| A = 2m=-oo J m (a)e lm>/ '). The above equations become: 
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We now define r F u , m (k z ,k ± ) = A J n dtf> i m F u (k)e ,m ^. This results in the following equations : 
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poo poo pm 

T u . m (k z ,k x ) = Idz rdr I d& E u (r) JJk ± r)e 
J — OO Jo Jo 

These equations correspond to Eqs. (6a) and (7a). 
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(A.6) 
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Appendix A.2. Cylindrical components 

Fet us now consider fields of the type E r , B, or J, , which we denote generally by F r . We have : 


F r = cos (9)F X + sin(0)F v = 


F x iF y w F x + iF y w 


(A.9) 


Using Eq. (A.8) leads to 
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where we relabeled the dummy variable m in the above sums. Let us thus define T . m = 1 + if y,m-i)/2 and 

T+sn = CF*,m+i - iT y , m+ 1 )/ 2 . This results in: 


| w s-*00 /~*00 

Fr(r) = ^2 2 J dk ~ J k ±dk ± (f'+.m Jm+\(k±r) + f-,m J m -\(k±rj)e 


imG+ik-z 


With the same definitions and the same method, it is also easy to show that: 

00 s-»oo s-»oo 

Fg(r ) = ^ Id k z I fc ± d/: ± i(f+, m J m +i(k±r) - l F_, m J m -\(k ± rf)e~" 
J -00 Jo 


im6+ik z z 
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Appendix B. Maxwell equations for the spectral coefficients 


In this section, let us derive the Maxwell equations for the spectral coefficients Eqs. (8a) to (8c) from the 
Maxwell equations written in cylindrical coordinates Eqs. (5a) to (5c). 

When replacing the Fourier-Hankel decomposition (Eqs. (6b), (6c) and (A.8)) in the Maxwell equations 
Eqs. (5a) to (5c), we first notice that the modes proportional to e~ lme+lkzZ for different values of in and k, are 
not coupled. These different modes can thus be treated separately. The same cannot be said of the modes cor¬ 
responding to different values of k ± , since they may be coupled through the Bessel functions J m (k±r ) and their 
derivatives. In the following, we write only the equations corresponding to d,B = -V x E, since the equation 
c~ 2 d,E — V x B - jjf)j can be treated very similarly. These equations become 


X oo 
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By taking the sum and difference of the first two equations, and by rearranging the third equation, we obtain: 
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We can now use the relations -r-.J m (k±r) + J'(k ± r) = J m -i(k ± r) and -P-J m (kj_r)-J'(k ± r) = J m+ i(k ± r ) (see relation 
9.1.27 in [51]), and obtain: 
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Each equation of the above system contains Bessel functions of only one given order (m + 1, m - 1 or m). This 
allows the different k ± components to be separated, since the functions J n (k ± r ), for a fixed n and different values 
of k ± , form a basis of the set of real functions: 


2 d t !B + , m = ik ± 8 z , m - 2 k z 8^ m 

(B.4a) 

2 — ik ± S zrn + lk z 8— m 

(B.4b) 

~ik±8+, m ik ± 8- m 

(B.4c) 


Appendix C. PSATD scheme, in the Fourier-Hankel representation 

We use a scheme very similar to that of [10]. In this scheme the currents are considered constant over one 
timestep, and the charge density is considered linear in time. 

Appendix C.l. Expressions for S m 

By combining Eqs. (8a) to (8c) and Eq. (9), one can find the propagation equations for B. 

d 2 S+, m + c 2 (k\ + k 2 )S+^ m = p 0 c 2 J z , m + k z J +itn j (C.la) 

d 2 A m + c 2 (k 2 ± + k 2 )A m = p 0 c 2 (- l ^J zjn - k z J ., m J (C.lb) 

d 2r A m + c 2 (k\ + k 2 ,)S Zitn = poc 2 (ik ± J+, m + ik ± J- tm ) (C. lc) 

Let us integrate these equations for f e [nAt, (n + l)Af], In this interval, J m (t) is constant and equal to A",, 
and thus the right-hand side of the above equations is constant. Using Green functions, the general solution of a 
differential equation of the form d 2 f + co 2 f = go, where go is a constant, is 

fit) = f(to) cos[ oXt - 1 0 ) ] + d t f(t o r [a)(t —— + ^(1 — cos[ oj(t — to) ]) (C.2) 

CO (jO 1 

We thus use the above expression, with co 2 = c 2 (k]_ + k ?), to integrate the fields from to = nAt to t - (n + l)Af. In 
particular, we use again the Maxwell equations Eqs. (8a) to (8c) to obtain the expression of d,!B m (to). This yields: 
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where C = cos(wAf) and S = sin(wAf). 
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Appendix C.2. Expressions for £>,„ 

Similarly, when combining Eqs. (8a) to (8c) and Eq. (9), the propagation equations for E are: 
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Let us again integrate these equations for t e [nAt, (n + 1)Ar]. In this interval, J m {t) is constant (thus its time 
derivatives drop), and p m is linear in time. As a consequence the right hand side is proportional to p n m + (p"f l - 
p" )(f ~ to) I At. Using Green functions, the solution of d 2 f + co 2 f = p n m + (p” +1 - p" n )(t - to)/At is 
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which, for t = to + At, reduces to 
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Appendix D. Discrete Hankel Transform 

Appendix D.l. Calculation of the transformation matrices M n m and M' nm 

Here, for the Discrete Hankel Transform, we use a transformation similar to that of [27, 25, 29], but we extend 
it to the case of an evenly-spaced grid in real space (as opposed to one that is distributed according to the zeros of 
the Bessel function, which would have been inconvenient for current deposition and field gathering). Moreover, 
we impose, as much as possible, that the succession of a Discrete Hankel Transform (DHT) and an Inverse Discrete 
Hankel Transform (IDHT) retrieves the initial function. As explained in the text of the article, we use a matrix 
formalism for the DHT: 


N r - 1 N r - 1 

DHT ”[/] (k'fj) - j] (M n , m )j' P f(r p ) IDHTH'Lg] {rf) - ^ {M' njn ) hp g(kf p ) (D.l) 

p =0 p =0 

where n is the order of the Hankel transform, and where m is the index of the spectral grid k”' . on which the 
Hankel transform is evaluated. In practice, as mentioned in the text, these transforms are only used in the cases 
n = m — 1, n = m or n = m + 1. 

The N r x N r matrices M n m and M' nm can be entirely determined by a set of N} constraints. These constraints 
can be found, for instance, by imposing the value of the DHT for a set of N r different functions, whose exact 
analytical Hankel transforms are known. In our case, we impose that the DHT be equal to the exact analytical 
Hankel transform for the eigenmodes of a cavity with perfectly conducting boundary at r max ( E(r max , z) = 0), since 
these physical eigenmodes should also be eigenmodes of our PIC cycle. These eigenmodes have the following 
form: 


E : 
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where © is the Heaviside function and where k'" f = a"'/r max , with a'" the fth positive zero of the Bessel function 
of order m. The exact Hankel transform of these modes, evaluated on the discrete set {k'f j) reads (see Appendix 
D.2 for a derivation) 


HT n [J n (k”f ( r)®(r max -r)](k^ j ) 


s2 , r; 

Jo 


rdr J n (kljr)J n (k^r) 


. [J„ +s ,Jaf)] 2 6j, t 
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where n is either m - 1, m or m + 1. Note that the above relation is valid for any value of m, j, f and n (provided 
that n e \m - 1, m, m + 1}), except when n 4- 0, m 4 0 and t - 0 simultaneously (again, see Appendix D.2 for an 
explanation). Here we impose that the DHT is consistent with Eq. (D.3), where applicable 


N r - 1 


DHT™ [ J n (k m u r)Q(r max - r) ] (kf f = ^ (M njn )j. p J n (k‘l { r p ) = ;r r 2 nax [J n+dnm (a™)] 2 5 jX (D.4) 


and we use the constraints given by Eq. (D.4) to obtain the matrix M n m 
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Case where n = 0 or m = 0. In this case, Eq. (D.4) is applicable for any j and £ in {0, ...,N r - 1), and thus this 
provides IV, 2 constraints on the matrix M nm , which allow one to completely determine it. In fact, a closer look at 
Eq. (D.4) shows that the inverse matrix M~ ] m can be directly extracted from the above relations, since this inverse 
matrices is defined by the relation Yip(^n,m)j,p(M~l n ) Pt ( = 6jx- Thus, from the above relations, one can directly 
infer: 


( M„jn) P ,c = 


Uk'l e r p ) 


nr max\- J n+6 n ., n (a n r ')] 


(D.5) 


The matrix M n m can then be extracted, by numerically inverting the matrices M" given by the above expression. 
In addition, we impose that M' m be exactly equal to M~ j rr so that the succession of a DHT and IDHT retrieves 
exactly the initial function. 


Case where n 4 0 and m 4 0. In this case, the equation Eq. (D.4) is not valid for t - 0, and thus it provides only 
N r (N r - 1) constraints on the matrix M„ m , which is not enough to completely determine it. In this case, we use an 
empirical method in which we impose the additional constraints 


(M n , m )o, P =0 for p e {0, ...,N r - 1) (D.6) 

This constraint is imposed because we choose the amplitude of the Hankel mode proportional to J n (k™r) to be 0 in 
the simulation. (For n 4 0 and m 4 0, J„(k' f "r) = 0 for any r, so that the amplitude of the corresponding mode has 
no physical meaning whatsoever.) This allows the matrix M njm to be entirely determined. In addition, to obtain 
Af' , we impose (M' nm ) p f = , " ±J p for any £ 6 - 1) (as in the case n = 0 or m = 0), but also 

(M' n m ) P '( = 0 for t = 0. This is done again for consistency with the fact that J n (k™r) — 0. We note that the above 
method gives satisfying results for m = 1 but not for m = 2. In the future, further work will be done towards a 
better Hankel transform representation. 


Appendix D.2. Derivation of Eq. (D.3) 

Let us first remark that, through a change of variable where r = r max t , Eq. (D.3) is equivalent to 

JJaffJfaJt) = [J n+s „Jcf”)] 2 8 jX (D.7) 

and let us prove this equation for different cases, excluding the case where n 4 0, m 4 0 and £ - 0 and in which it 
is not valid. 



Case where n = m and £4 0. In this case, we use the relation 11.4.5 of [51]. Since n = m, and since by definition 
a'J is the zero of the Bessel function of order m, we have J n {a"j) = J„(a" r ) = 0, and the relation 11.4.5 is thus used 
in the case b = 0 and a - 1 (see relation 11.4.5 in [51] for the definition of the a and b coefficients). This yields 

ftdt J„(ap)J„(aJt) = l -U' n (a™)] 2 d j , e = i[7„ + i(a”)] 2 ^ (D.8) 

where we further used the relation 9.1.27 in [51], which reads J' n (a"') = -J n+ \{a"‘) + and took into 

account the fact that J n (a"') = 0 in our case. 


Case where n 4 m and £4 0. In this case, n is either m + 1 or m — 1, since we restricted ourselves to n e 
{in - l,m,m +1). Let us prove the equation Eq. (D.7) in the case n - m + 1 (the proof for n — m - 1 being very 
similar). By definition, we have J nt (a n ( ') = 0, and thus from the relation 9.1.27 in [51] 7' M+1 (a'")+ ^r7 m+ i(a“) = 0. 
We can thus apply the relation 11.4.5 in [51] with a = m + 1 and b = 1, and find 


f tdt JfapVniaJt) = —4a { (m + D 2 + Kf - ” 2 ) = hjnWfSjt 

Jo ^ 


(D.9) 


Case where m = 0 and £ = 0. The proof of the two above cases used the relation 11.4.5 from [51], which is 
applicable as long as a'" > 0. This is the case for £40 (i.e. the two above cases) but also for m = 0 and £ — 0 (i.e. 
Jo(0) 4 0 and thus a g >0). Asa consequence, the above proofs are also valid for m = 0 and £ = 0. 
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Case where n = 0, m + 0 and (' = 0. In this case, a'" = 0 and thus the relation 11.4.5 from [51] does not apply. 
However, for j + 0, the left-hand side of equation Eq. (D.7) reduces to 


Jtdt J n (a™t)J n (a n ‘t) = J 0 ( 0) J tdtJ n (a'Jt ) = ^ J )tdtJ n (t) = 0 (D.10) 


where we used relation 11.1.1 in [51]. On the other hand, for j — 0, one has a 1 " = a'." = 0 and thus 


Jtdt J n (ap)J n (a n ;t) = [7 0 (0)] 2 Jtdt = |[-/o(0)] 2 


(D.l 1) 
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